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ABSTRACT 

We present the results of radiation-magnetohydrodynamic simulations of the formation and 
expansion of H II regions and their surrounding photodissociation regions (PDRs) in turbulent, 
magnetised, molecular clouds on scales of up to 4 parsecs. We include the effects of ionising 
and non-ionising ultraviolet radiation and x rays from population synthesis models of young star 
clusters. For all our simulations we find that the H II region expansion reduces the disordered 
component of the magnetic field, imposing a large-scale order on the field around its border, 
with the field in the neutral gas tending to lie along the ionisation front, while the field in the 
ionised gas tends to be perpendicular to the front. The highest pressure compressed neutral 
and molecular gas is driven towards approximate equipartition between thermal, magnetic, and 
turbulent energy densities, whereas lower pressure neutral/molecular gas bifurcates into, on 
the one hand, quiescent, magnetically dominated regions, and, on the other hand, turbulent, 
demagnetised regions. The ionised gas shows approximate equipartition between thermal 
and turbulent energy densities, but with magnetic energy densities that are 1 to 3 orders of 
magnitude lower. A high velocity dispersion (~ 8 km s" 1 ) is maintained in the ionised gas 
throughout our simulations, despite the mean expansion velocity being significantly lower. 
The magnetic field does not significantly brake the large-scale Hll region expansion on 
the length and timescales accessible to our simulations, but it does tend to suppress the 
smallest-scale fragmentation and radiation-driven implosion of neutral/molecular gas that 
forms globules and pillars at the edge of the H II region. However, the relative luminosity of 
ionising and non-ionising radiation has a much larger influence than the presence or absence 
of the magnetic field. When the star cluster radiation field is relatively soft (as in the case of 
a lower mass cluster, containing an earliest spectral type of B0. 5), then fragmentation is less 
vigorous and a thick, relatively smooth PDR forms. Accompanying movies are available at 
http : //youtube . com/user/divBequalsO 

Key words: H II regions - ISM: kinematics and dynamics - magnetohydrodynamics - photon- 
dominated region (PDR) - Stars: formation 



1 INTRODUCTION 

H II regions, that is, regions of photoionised gas, are among the most 
arresting astronomical objects at optical wavelengths. The basic the- 
ory behind their formation and expansion has been known for some 
time (Stromgren 1939; Kahn 1954; Spitzer 1978). More recently, 
attention has turned to explaining the irregular structures, filaments, 
globules and clumps seen within and around these regions. These 
could be due to underlying density inhomogeneities in the molecular 
cloud into which the H II region is expanding, or could be due to 
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instabilities at the ionisation front itself (Garcia-Segura & Franco 
1996; Mellema et al. 2006a; Williams et al. 2001; Gritschneder et 
al. 2010; Mackey & Lim 2010a). The effect of the stellar ionising 
radiation on these structures can lead to radiation driven implosion, 
and the subsequent compression could result in the formation of 
new stars (Bertoldi 1989; Esquivel & Raga 2007; Motoyama et al. 
2007). 

Magnetic fields are well known to pervade molecular clouds 
and are thought to play an important role in regulating the star- 
formation process by providing support to collapsing structures in 
molecular clouds (Mouschovias & Ciolek 1999), although recently 
their relative importance in this process has become less clear (Mc- 
Kee & Ostriker 2007). However, theoretical work on the interplay 
between ionising radiation and magnetic fields in and around H II 
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regions has only been possible recently, with the advent of radiation- 
magnetohydrodynamic codes. Krumholz et al. (2007) performed 
the first simulations of the expansion of an H II region in a uniform, 
magnetised medium. At early times, the thermal pressure of the pho- 
toionised gas dominates over the magnetic field, which is pushed 
out of the expanding H II region. At late times (several Myr), the 
photoionised region becomes elongated along the field lines and 
the magnetic field filled back in. Henney et al. (2009) and Mackey 
& Lim (2010b) have studied the effects of ionising radiation on 
magnetised globules, finding that the photoevaporation process is 
significantly altered only when the magnetic pressure exceeds the 
thermal pressure by more than a factor of ten in the initial globule. 
Peters et al. (2010) have recently studied the combined effects of 
magnetic fields and photoionisation on the formation of a massive 
star on scales < 0. 1 pc, finding that rotational winding of the field 
lines produces a magnetised 'bubble', whose magnetic pressure acts 
to confine the nascent H II region during its ultracompact phase. 

In this paper, we study the expansion of H II regions in non- 
uniform (turbulent) magnetised molecular clouds. Previously, we 
studied H II region expansion in a turbulent medium without a mag- 
netic field and our results were strikingly similar to observed H II 
regions (Mellema et al. 2006a). In the present work, we study the 
effect of the expanding photoionised gas on structures in the sur- 
rounding magnetised medium and also the effect that the turbulent 
magnetic field has on the H II region itself. At least at the relatively 
early times studied in our simulations, most of the interesting ef- 
fects due to the magnetic field will occur within the neutral, photon 
dominated region around the H II region, and so we use a careful 
treatment of the heating and cooling processes in this gas, first dis- 
cussed in Henney et al. (2009), in order to have a realistic portrayal 
of the dynamical processes here. 

The wealth of new observational data at infrared wavelengths 
enables comparisons to be made of our simulations not only with 
observations of the photoionised gas, as done previously in Mellema 
et al. (2006a) but also with the PDR and molecular gas affected by 
the H II region. In particular, the Spitzer GLIMPSE surveys of bub- 
bles in the Galactic plane (Churchwell et al. 2006, 2007) highlight 
polycyclic aromatic hydrocarbon (PAH) emission at the outer edges 
of H II regions and in PDRs. These surveys show that the thickness 
of the PAH-emitting shell varies between 0.2 and 0.4 of the outer 
radius and that shell thickness increases approximately linearly with 
radius. Moreover, the bubbles are generally non- spherical. Herschel 
and APEX observations reveal thermal emission from cold dust 
in the molecular gas outside the PDR (e.g. Anderson et al. 2010; 
Deharveng et al. 2010). These observations are useful in identifying 
condensations that could be collapsing to form new stars. 

Studies of magnetic fields in and around molecular clouds are 
important for understanding the star- formation process. There is 
no consensus as to the importance of magnetic fields in the forma- 
tion of clouds, cores, and ultimately stars. If magnetic fields are 
dynamically important, then star formation is regulated by ambipo- 
lar diffusion (e.g., Mouschovias & Ciolek 1999). Alternatively, if 
magnetic fields are unimportant then other processes such as turbu- 
lence and stellar feedback provide support to molecular clouds and 
determine the star-formation efficiency (e.g., Ballesteros-Paredes 
et al. 1999; Elmegreen 2000; Matzner 2002; Krumholz, Matzner, 
& McKee 2006; Tasker & Tan 2009; Tasker 201 1). It is not known 
whether magnetic fields are strong enough to dynamically support 
molecular clouds and more evidence is needed. It is possible that 
reality lies somewhere between these extremes. Recent Zeeman 
observations of molecular clouds (Crutcher, Hakobian, & Troland 
2009; Crutcher et al. 2010) do not agree with predictions of ideal- 



ized ambipolar diffusion models (strong magnetic field) but more 
observations are needed. 

There are three main techniques for tracing the magnetic field 
in diffuse HI and molecular clouds: polarisation from aligned dust 
grains, linear polarisation of spectral lines, and Zeeman splitting of 
spectral lines (Heiles & Crutcher 2005). Polarisation studies yield 
B ± , the strength of the field projected onto the plane of the sky, 
whereas Zeeman splitting gives B\\, the strength of the field parallel 
to the line of sight. In particular, the H I Zeeman effect has been 
used to study the atomic gas (assumed to be in the photodissoci- 
ated region) in the star-forming regions Ml 7 (Brogan & Troland 
2001) and the Ophiuchus dark cloud complex (Goodman & Heiles 
1994), while the OH Zeeman effect is used to study the line-of-sight 
magnetic field in the molecular gas. The results of such studies 
suggest that there is rough equipartition of magnetic and turbulent 
(dynamical) energy densities in the neutral gas. In the Ophiuchus 
dark cloud region, the uniform field in the H I gas is estimated to be 
10.2 iiG (Goodman & Heiles 1994), whereas the magnitude of the 
line-of-sight magnetic field strength in the H I gas associated with 
M17 is found to be ~ 500//G (at 60" resolution, or even twice as 
high at 26" resolution) (Brogan & Troland 2001), suggesting that 
there is small-scale structure in the magnetic field. Recent linear 
polarisation studies at submillimetre wavelengths of the magnetic 
field around the ultracompact Hll region G5. 89-0.39 (Tang et al. 
2009) provide evidence for compression of the B field in the sur- 
rounding molecular cloud and general disturbance of the magnetic 
field caused by the expansion of the H II region. 

For photoionised gas, Faraday rotation measures of line-of- 
sight extragalactic radio sources are used to calculate the magnetic 
field strengths in foreground H II regions (Heiles et al. 1981). The 
Zeeman effect cannot be used in the ionised gas because the thermal 
linewidths are too large. Heiles et al. (1981) use double extragalactic 
radio sources to sample two positions in the Hll regions SI 17 and 
S264 and find identical rotation measures in each case of order 1 
or 2 nG. In these regions, the thermal energy density dominates the 
magnetic energy density. In a third H II region, SI 19, for which only 
a single rotation measure is available, the line-of-sight magnetic 
field strength is around 20 //G and the ratio of magnetic to thermal 
energy density is around 0.4, hence the magnetic field could have 
had an effect on the expansion of this H II region. 

The remainder of this paper is organised as follows. In § 2 we 
describe our numerical algorithm and present two test problems: 
the classical Spitzer law for the expansion of a non-magnetised H II 
region in a uniform medium, and the expansion of an H II region in a 
uniform magnetised medium (Krumholz et al. 2007; Mackey & Lim 
2010b). In the second case, we draw attention to some interesting 
properties of the solution that have not been described previously. 

In § 3 we describe the expansion of H II regions in both mag- 
netised and non-magnetised turbulent media, taking as initial con- 
ditions for the ambient medium the results of an MHD turbulence 
calculation (Vazquez-Semadeni et al. 2005), and considering both 
a strong (approximately 09 spectral type) ionising source and a 
weak (B0.5) one. In § 4 we make qualitative comparisons between 
simulated optical and long- wavelength images and observations. We 
also produce synthetic maps of the projected line-of-sight and plane- 
of-sky components of the magnetic field, which can be qualitatively 
compared with observations of the magnetic field in star-forming 
regions. We discuss the globules and filaments at the periphery of 
our simulated H II regions in the context of recent numerical work 
on photoionised magnetised globules (Henney et al. 2009; Mackey 
& Lim 2010b). Finally, in § 5 we outline our conclusions. 



© 0000 RAS, MNRAS 000, 000-000 



Radiation-MHD simulations ofHil regions and PDRs 3 



2 NUMERICAL MODEL 



2.1 Numerical Method 

The numerical method used in this paper is the same as that de- 
scribed by Henney et al. (2009). The Phab-C 2 code combines the 
ideal magnetohydrodynamic code described by De Colle & Raga 
(2006) with the C 2 -Ray method (Conservative Causal Ray) for ra- 
diative transfer developed by Mellema et al. (2006b), and adds a 
new treatment of the heating and cooling in the neutral gas (Henney 
et al. 2009). The MHD code is a second-order upwind scheme that 
integrates the ideal MHD equations using a Godunov method with a 
Riemann solver, similar to that described by Falle et al. (1998) but 
with the constrained transport (CT) method (e.g., Toth 2000) used to 
maintain the divergence of the magnetic field close to zero. An artifi- 
cial (Lapidus) viscosity is also included (Colella & Woodward 1984) 
in order to broaden the discontinuities and reduce oscillations. The 
C 2 -Ray method for radiative transfer is explicitly photon conserving 
(Mellema et al. 2006b) and the algorithm allows for timesteps much 
longer than the characteristic ionisation timescales or the timescale 
for the ionisation front to cross a numerical grid cell, by means of 
an analytical relaxation solution for the ionisation rate equations. 
This results in a highly efficient radiative transfer method. 

The MHD and radiative transfer codes are coupled via operator 
splitting, the inclusion of equations for the advection of ionised and 
neutral density in the MHD code, and the energy equation, which 
includes both radiative cooling and heating due to the absorption of 
stellar radiation by gas or dust: ionising extreme ultraviolet radiation 
in the photoionised region, non-ionising far ultraviolet radiation in 
the neutral gas and x-rays in the dense molecular gas. A detailed dis- 
cussion of the heating and cooling contributions is given by Henney 
et al. (2009). In particular, the treatment of heating in the neutral 
PDR region, taking into account FUV/optical dust heating and x- 
ray heating, is a substantial improvement over standard treatments 
in the literature. The heating term for the neutral gas is calibrated 
using Cloudy (Ferland et al. 1998) and is tailored to the ionisation 
source. Moreover, x-rays and FUV radiation will be enhanced due 
to the presence of associated low-mass stars in the star-formation 
environment, and we use the results of Fatuzzo & Adams (2008) 
to estimate the distributions of FUV and EUV fluxes appropriate 
to stellar clusters of membership sizes relevant to the two cases we 
study: an 09 and a B0. 5 star. The energy equation is solved itera- 
tively using substepping and the timestep for the whole calculation 
is determined solely by the magnetohydrodynamic timestep. 

The individual components of the Phab-C 2 code have been 
tested against standard problems (De Colle 2005; De Colle & Raga 
2004; Iliev et al. 2006) and the whole code was applied to the prob- 
lem of the photoionisation of magnetised globules Henney et al. 
(2009), work which has recently been independently corroborated 
by Mackey & Lim (2010b). In this paper, we have implemented 
an entropy fix (Balsara & Spicer 1999) to correctly deal with the 
thermal and dynamical pressures in regions of low density where the 
magnetic pressure dominates. We have also improved the treatment 
of the boundary conditions because they can become important, 
particularly towards the end of the simulations. For our O star simu- 
lations we use outflow boundary conditions, since thermal pressure 
is expected to dominate throughout the evolution and the H II region 
will rapidly grow to the size of the computational box. For the B 
star simulations, the expansion timescale of the H II region is much 
longer and is comparable to the crossing time of the computational 
box corresponding to the initial turbulent velocity dispersion of the 
molecular gas. We therefore use periodic boundary conditions in 
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Figure 1. Radius against time for the hydrodynamic expansion of a pho- 
toionised region in a uniform medium. Bottom panel: numerical simulation 
(black line) and analytical solution (gray line). Top panel: relative error. 

this case, which ensures that the total mass in the computational box 
remains constant. 

Parallelisation using Open-MP was implemented where possi- 
ble, enabling the current simulations to be run in about two weeks 
on an Intel 8 -core server with 32 GB of RAM. Preliminary work 
was performed on the Kan Balam supercomputer of the Universidad 
Nacional Autonoma de Mexico. 

2.2 Test Problems 

2.2. 1 Hydrodynamic Hll Region Expansion in a Uniform Medium 
Org/ 

The radiation part of the code has been extensively tested and 
documented (Mellema et al. 2006b; Iliev et al. 2006). In order to 
test the radiation (magneto-)hydrodynamics combined code we first 
consider the purely hydrodynamic evolution of a photoionised region 
around an ionising source in a uniform medium. The analytical 
solution to this problem is the familiar Spitzer (1978) expansion law 

^Spitzer =tf |l + ~\ , (1) 

where /^spitzer is the ionisation front radius, Rq = (3 (2^47™;^) 1/3 is 
the initial Stromgren radius, c x is the sound speed in the ionised gas 
and t is time, <2h is the ionising photon rate, n is the ambient density 
and a B is the case B recombination rate, a B = 2.6 x io~ 13 r -0 - 7 (Os- 
terbrock & Ferland 2006), where T is the temperature in the ionised 
gas. For an ambient density n = 1000 cm -3 , stellar ionising photon 
rate Q n = 5 x 10 46 s" 1 and stellar temperature T cS = 37500 K the 
results are shown in Fig. 1, together with the corresponding analytic 
solution. These parameters were chosen because they represent the 
mean values of the turbulent medium B star case described below 
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Figure 2. Mean gas velocity against time for the hydrodynamic expansion 
of a photoionised region in a uniform medium. Black line — numerical 
simulation; thick gray line — analytic solution. 

in § 3. The lower panel shows radius against time, while the upper 
panel shows the relative error, defined as 

_, ^ion — ^Spitzer / _ x 

Error = . (2) 

^Spitzer 

From this figure, we see that the largest discrepancies, of about 
6%, occur at the beginning of the expansion, when a shock wave 
begins to propagate outwards ahead of the ionisation front and a 
corresponding rarefaction wave travels back towards the ionising 
source, which is not accounted for in the analytical solution. At 
later times, the error is never more than 1%. Our code, therefore, 
compares extremely favorably with other radiation-hydrodynamics 
codes on this problem (see, e.g. Arthur & Hoare 2006; Krumholz 
et al. 2007; Iliev et al. 2009). There is no evidence for overcooling 
in the ionisation front in our simulations, as was found to be a 
problem by Krumholz et al. (2007). In the analytic Spitzer solution 
we assume an H II region temperature T = 8900 K, which gives 
the best fit to our numerical result. This is about 400 K hotter than 
the volume-averaged ionised gas temperature in our simulation, 
but is roughly equal to the ionised gas temperature just inside the 
ionisation front. 

In Fig. 2 we plot the gas radial velocity from the numerical 
simulation and the corresponding analytical solution. The numerical 
solution closely follows the analytical solution, though with oscilla- 
tions due to acoustic phenomena, such as the initial rarefaction wave, 
which propagates back into the H II region when the ionization front 
begins to expand outwards. 

2.2.2 MHD H II Region Expansion 

There is no corresponding analytical solution for the expansion of 
an H II region into a magnetised medium. Therefore, in order to test 
the MHD part of the numerical code, we carry out a simulation of 
the expansion of an H II region in a uniform magnetised medium 
using the same parameters as those of Krumholz et al. (2007). This 
problem was also used to test the radiation MHD code of Mackey & 



Lim (2010b). The simulation is performed in a 25 6 3 computational 
box with constant ambient density and temperature of n = 10 2 cm -3 
and T = 1 1 K, respectively, a constant ambient magnetic field 
directed at 30° to the x-axis in the x-y plane 1 of strength B = 
14.2//G and an ionising source producing <2h = 4 x 10 46 ionising 
photons s _1 . The magnetic field is placed at an angle to the grid lines 
in order to be able to distinguish between physical and computational 
effects. We run the simulation in a (20 pc) 3 box and in a (40 pc) 3 
box, the former to study the early evolution (t < 2 Myr) and the 
latter to highlight the late-time evolution. We note that the ionising 
source is very weak (later than spectral type B0.5, Vacca et al. 1996) 
and that the ambient density is lower than that typically found in 
star-forming regions. 

Krumholz et al. (2007) found a critical radius and time for 
when magnetic pressure and tension start to become important for 
the expansion of the H II region in a uniform magnetised medium, 
which corresponds to when the magnetic pressure in the neutral gas 
becomes of the same order as the thermal pressure in the ionised 
gas, i.e., p n v 2 A ~ p x c\. Here, p n , p x are the densities in the neutral 
and ionised gas, v A is the Alfven speed, and c x is the sound speed in 
the ionised gas. Thus, magnetic effects become significant when the 
H II region has expanded to a radius of roughly 

R m = l-M R , (3) 

where R is the initial Stromgren radius in a non-magnetised medium 
of the same uniform density. The corresponding time is 

4/ Ci \ 7/3 

t^j(-) k. (4) 

where t is the sound crossing time of the initial Stromgren sphere, 
and both of these equations implicitly assume that the classical 
Spitzer expansion law (Spitzer 1978) holds until the magnetic field 
becomes important. 

For the parameters in this test problem, the initial Stromgren 
radius is R = 0.48 pc, and the Alfven speed in the ambient magne- 
tised medium is 2.6 km s -1 . Our simulations give a temperature in 
the photoionised gas of T x ~ 9000 K and hence a sound speed of 
ci = 9.8 km s _1 (compare the lower assumed ionised sound speed 
c\ ~ 8.7 km s _1 in the Krumholz et al. 2007 models). Thus, the criti- 
cal radius and time are R m = 2.8 pc and t m = 0.61 Myr, respectively, 
for our models. 

In Figs. 3 and 4 we show snapshots of the magnetic field and 
the photoionised gas distribution in the central x-y, x-z and y-z 
planes of the computational box after 2 and 6 Myr of evolution. 
The coloured discs shown between the panels of the top row are 
keys to the interpretation of the magnetic field images. The hue 
indicates the projected field orientation in the plane of each cut (red 
for vertical, cyan for horizontal, etc.), while the brightness indicates 
the magnetic field strength, as is shown in the left-hand disc. The 
colour saturation diminishes as the out-of -plane component of the 
field increases, as shown in the right-hand disk. This non-standard 
representation allows the magnetic field to be sampled at every pixel. 
The lower panels of the figures show the ionisation fraction, temper- 
ature and density of the gas in the same central planes. Ionisation 
fraction is represented by colour: red indicates fully ionised gas 
while blue/green is partially ionised. The gas temperature varies be- 
ween ~ 10 4 K (red) in the ionised gas to a few hundred Kelvin in the 

1 The original simulation performed by Krumholz et al. (2007) has the 
magnetic field parallel to the x-axis. 
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Figure 3. Expansion of an H II region in a magnetised uniform medium after 2 Myr in a (20 pc) 3 computational box. The top row shows the magnetic field in 
cuts in the central x-y (left), x-z (centre) and y-z (right) planes. The colour indicates the direction of the magnetic field, where the coloured discs show the field 
angle in the plane for a strong field (left, brighter disc) and a weaker field (right, dimmer disc). Grey indicates magnetic field perpendicular to the plane. The 
bottom row shows ionisation fraction x, temperature T and density p, coded as 'hue', 'saturation', and 'value', respectively. Red indicates hot (~ 10 4 K) ionised 
gas, blue/green is partially ionised gas, purple is warm (~ 300 K) neutral (PDR) gas and grey indicates cold neutral (molecular) gas. The density is represented 
by the intensity, with the densest gas being white and the most diffuse is black. 



neutral (PDR) gas. The gas density is represented by the intensity in 
these figures, with the densest, cold molecular gas appearing white. 

At early times, t < t m , the Hll region is essentially spherical 
and expands at approximately the same speed (half the sound speed) 
in all directions because the ionised gas thermal pressure is much 
higher than the ambient magnetic pressure. Expansion across the 
magnetic field lines compresses the B-field. The expanding Hll 
region and PDR are preceded by a fast-mode shock, which travels 
furthest perpendicular to the magnetic field direction. The fast-mode 
shock bends and compresses the field lines and so the PDR and 
H II region propagate in a modified magnetic field. The slow-mode 
shock is strongest along the field lines and most evident at the end 
caps, where it produces a dense shell. It compresses the gas and 
demagnetises it (reduces the field). The slow-mode shock is almost 
isothermal and the dense shell behind it is pushed by the pressure in 
the photoionised region. At early times the magnetic field is reduced 
in the photoionised gas by the expansion of the H II region. At later 
times, t > t m , the magnetic tension pulls the field back into the Hll 
region at the equator. 

In Fig. 3 we show the situation after 2 Myr, when the magnetic 
field has started to pull back into the H II region. The H II region 
is elongated along the magnetic field direction but the photon dom- 



inated region (warm neutral gas) is approximately spherical. The 
magnetic field retains its original direction except around the edge 
of the H II region. This is particularly evident at the poles of the H II 
region. At the time shown in the figure, the magnetic field appears 
concentrated towards the centre, with regions of lower magnetic 
field around the equator or waist of the structure. Neutral gas is 
pulled into the H II region along with the magnetic field, becomes 
ionised and flows away from the ionisation front into the H II region. 
It then flows along the field lines and is channelled out towards the 
poles, where the ionisation front transforms into a recombination 
front. We note that at this time in the simulation, we do see insta- 
bilities form in the direction perpendicular to the magnetic field, as 
reported by Krumholz et al. (2007) and Mackey & Lim (2010b). 
These instabilities form at the equator where the ionisation front is 
parallel to the magnetic field direction because the slow-mode shock 
which precedes the ionisation front, is not permitted to cross the 
field lines. As a result, the ionisation front corrugates. 

Our Fig. 4 shows the simulation after 6 Myr in a (40 pc) 3 com- 
putational box. The H II region has stopped expanding perpendicular 
to the magnetic field but the shocked neutral shell is still expanding 
along the field lines, leading to a very elongated structure. However, 
in the direction along the field lines, the flow has become one dimen- 



© 0000 RAS, MNRAS 000, 000-000 



6 S. J. Arthur et al. 





i 




o 






I* 


t = 6.00 Myr 



Figure 4. Same as Fig. 3 but after 6 Myr in a (40 pc) 3 computational box. 



sional and so the density does not fall off as the H II region expands. 
Thus, there is a limit as to how far the ionisation/recombination front 
can travel in this direction. There is an extended region of partially 
ionised gas beyond the recombination front. The warm, neutral gas 
(PDR) is still roughly spherical in shape but is pierced by the cigar- 
shaped H II region along the field direction. The dense neutral knots 
formed as a result of the instability at the equator shadow parts of 
the PDR. The knots formed closest to the perpendicular direction 
remain there throughout the simulation, pinned by the magnetic 
field. The knots slightly off the equator experience the rocket effect 
due to the photoionisation of their outer skins. They move like beads 
on a wire along the magnetic field lines around the edge of the H II 
region, a process which can be clearly observed in animations of 
this simulation. By the end of the calculation, the magnetic field has 
refilled the H II region and become roughly uniform and returned 
completely to its original direction. 

Finally, in Figs. 5 and 6, we examine the expansion of the 
magnetised H II region. Fig. 5 shows the evolution of the radius of 
the photoionised region with time. The mean radius is consistently 
below the classical Spitzer expansion law, unlike the non-magnetic 
case (see Fig. 1). In the direction perpendicular to the magnetic 
field, the expansion slows with time, as shown by the evolution 
of the minimum radius. In the direction along the field lines, the 
expansion initially (t < 1 Myr) follows the classical Spitzer law but 
then accelerates as material is channelled along the field lines and 
out of the poles of the H II region. This expansion abruptly stops 
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Figure 5. Radius against time for the magnetohydrodynamic expansion of 
a photoionised region in a uniform magnetised medium. Numerical MHD 
simulation (black line) and analytical non-magnetised solution for classical 
Spitzer H II region expansion (gray line). Also shown are the polar radius 
(dotted line) and equatorial radius (dashed line). The break in the lines at 
t = 2 Myr is a result of plotting velocities from two simulations at different 
numerical resolutions. 



just after 4 Myr, which corresponds to the time when the ionisation 
front can no longer expand along the field lines because the flow has 
become one dimensional. 
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Figure 6. Mean velocities against time for the magnetohydrodynamic ex- 
pansion of a photoionised region in a uniform magnetised medium. Rms 
velocity (black line) and analytical non-magnetised solution for classical 
Spitzer Hll region expansion (thick grey line). Also shown are the mean 
radial velocity in the ionised gas (dashed line) and the velocity of the gas 
at the ionization front (dotted line). The break in the lines att = 2 Myr is 
a result of plotting velocities from two simulations at different numerical 
resolutions. 



In Fig. 6 the mean radial and root-mean-square (rms) velocities 
of the ionised gas are plotted together with the analytical value and 
also the mean velocity of the ionisation front. The discontinuity 
at 2 Myr is due to the change in resolution of the simulation. The 
calculated values are all above the analytical values and indicate 
that typical velocities in the photoionised gas are of the order 1- 
2 km s _1 throughout the simulated lifetime of the Hll region. The 
mean radial velocity shows the same sort of ringing that we saw in 
the non-magnetised simulation (see Fig. 2). After about 3 Myr, the 
gas velocities become dominated by the one-dimensional flow being 
channelled along the field lines. 



3 H II REGION EXPANSION IN A MAGNETISED 
TURBULENT MEDIUM 

In this section we describe the results of radiation MHD simulations 
of H II region expansion in a more realistic medium where neither 
the initial density distribution nor magnetic field are uniform. We 
consider both an O star ionising flux and one more appropriate to a 
B star in order to study the effect of the magnetic field on the ionised 
gas and on the neutral gas of the PDR. 



3.1 Initial Set Up 

As initial conditions we take the results of a 25 6 3 driven-turbulence, 
ideal magnetohydrodynamic simulation by Vazquez-Semadeni et 
al. (2005), specifically the simulation with initial isothermal Mach 
number M s = 10, Jeans number J - A and plasma beta p - 0.1. 
The original Vazquez-Semadeni et al. (2005) simulations are scale 
free and are characterised by the three nondimensional quantities 
M s = cr/c (the rms sonic Mach number of the turbulent velocity 
dispersion cr), / = L/Lj (the Jeans number, giving the size of the 
computational box L in units of the Jeans length Lj), and the plasma 
beta,/? = 8^p c 2 /5q (ratio of thermal to magnetic pressures). 

We set the scaling by choosing the computational box size to 
be 4 pc, giving the ambient temperature as T ~ 5 K, the mean 
atomic number density as (hq) = 1000 cm -3 and the mean magnetic 



field strength as B = 14.2yi/G. We choose as our starting point the 
time in the evolution when there is one collapsing object. At this 
point, the mean values of the magnetic field the number density 
are unchanged, since flux and mass are conserved, while the mean 
plasma beta is now p = 0.032, which is consistent with the rms 
value of the magnetic field = 24.16yt/G. That is, the turbulence 
has enhanced the rms magnetic field by a factor rj = 24.16/14.2, 
which leads to a change by a factor rf ~ 3 in the value of p over the 
original, uniform conditions. 

As in our previous paper (Mellema et al. 2006a), we place our 
ionising source in the centre of the collapsing object, remove 32M© 
of material (corresponding to the mass of the star formed) and take 
advantage of the periodic boundary conditions of the turbulence sim- 
ulation to move the source to the centre of the grid. We also subtract 
the source velocity from the whole grid. We consider two ionising 
sources: one corresponds approximately to an 09 star (Vacca et 
al. 1996), having an effective temperature r eff = 37500 K and an 
ionising photon rate of g H = 10 48 5 photons s _1 , while the other 
ionising source has <2h = 5 X 10 46 photons s _1 , which corresponds 
to B0.5. 



3.2 Results 

In this section we present our results both for the O star and the 
B star and in each case for the magnetic and also for the purely 
hydrodynamic simulations. In this way we can assess the importance 
of the magnetic field and the ionising photon flux on the evolving 
structures in the expanding Hll regions and surrounding neutral gas. 
We first present emission-line images of the simulated H II regions 
for both O and B stars in the magnetic and pure hydrodynamic 
cases when they have reached a similar size, which corresponds to 
200,000 yrs of evolution for the O star case and 10 6 yrs of evolution 
in the B star case. This permits a direct comparison with observations 
at optical and longer wavelengths. We then examine the expansion 
and time evolution of the ionised, neutral and molecular components 
using global statistical properties. The importance of the magnetic 
field in the different components is then studied for the particular O 
and B star simulations presented earlier. 

Although the ionised/neutral transition of hydrogen is treated 
explicitly in a detailed manner in our radiation-MHD code, the 
neutral/molecular transition is treated much more approximately. 
We simply assume that the molecular fraction is a predetermined 
function of the extinction A v from the central star to each point in 

our simulation: n mo \/n neut = (l + e~ A{Ay ~ 3) ^j , which was determined 
by fitting to the Cloudy models described in the Appendix of Henney 
et al. (2009). 

A rough estimate of the importance of the magnetic field in 
our turbulent media simulations can be obtained by calculating the 
critical radius and timescale, R m and t m (see Eqs. 3 and 4), for 
the initial value of the magnetic field and the mean density. Using 
#rms - 24.16 jiG and (n) = 10 3 cm -3 , and assuming a mean particle 
mass of 1.3 m H , the initial representative Alfven speed in the neutral 
gas is v A = 5 rms /(4^-p) 1/2 1.46 km s _1 . In the photoionised gas, 
the sound speed is c x = 9.8 km s _1 . For the strong ionising source 
(O star), the Stromgren radius in a uniform medium of density equal 
to the mean density would be Ro = 0.45 pc, giving a critical radius 
R m = 5.7 pc and time t m = 2.2 Myr for when the magnetic field 
starts to become important. If, instead of the rms magnetic field we 
use the mean magnetic field, (B) = 14.2//G, then the corresponding 
values are R m = 1 1.86 pc and t m = 7.9 Myr. For the weak ionising 
source (B star), the numbers are R = 0.113 pc, R m = I A3 pc 
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(R m = 2.98 pc) and t m = 0.55 Myr (t m = 1.98 Myr) depending on 
whether the rms or mean magnetic field is considered. 

Thus, for the strong ionising source, the critical diameter is 
larger than the dimension of the computational box, hence we do 
not expect that the magnetic field would have much effect on the 
expansion of the H II region in this case. However, for the weak 
ionising source, it is not clear whether the rms or the mean value 
of the magnetic field is most important for determining the global 
evolution of the Hll region. If the rms field is the determining 
parameter then the fact that the critical diameter in this case is 
smaller than the size of the computational box leads us to expect 
that the magnetic field will affect the global expansion of the H II 
region. On the other hand, if it turns out that the mean magnetic 
field is most important, then the effect will not be noticed within 
our computational size and timescales. Furthermore, given the non- 
uniform initial conditions, the timescale for the magnetic field to 
become important could vary depending on the local gas sound 
speed to Alfven speed ratio. 

3.2.1 Morphologies and Images 

We start our comparison by considering the morphological appear- 
ance of the H II regions and the surrounding neutral and molecular 
material. In Figs. 7 and 8 we show two types of image data, one 
showing mostly the ionised gas, where the different colours repre- 
sent emission from optical emission lines, the other shows synthetic, 
long- wavelength (infrared and radio) emission, emphasizing the 
neutral and molecular material. 

The optical emission was calculated as described in our previ- 
ous papers (Mellema et al. 2006a; Henney et al. 2009), assuming 
heavy element ionisation fractions that are fixed functions of the 
hydrogen ionisation fractions. For the B star models, these functions 
were recalibrated for the softer stellar ionising spectrum using the 
Cloudy plasma code (Ferland et al. 1998). For the O star models, 
we employ the "classical" HST red-green-blue colour scheme of 
[Nil], Ha, [O III]. Since the [O III] emission from the B star models 
is predicted to be very weak, in this case we use the scheme [S II], 
[Nil], Ha. In both schemes, the progression from red through green 
to blue corresponds to increasing degree of ionisation inside the H II 
region. The emission from all the optical lines is negligible in the 
dense neutral/molecular zones, but the dust absorption associated 
with these dense regions is visible in the images. Scattering by dust 
is not included. 

For the long-wavelength images the emission bands were cho- 
sen to give a more global view of the simulations than can be pro- 
vided by the optical emission lines. The red band shows the total 
column density of neutral/molecular gas, which very crudely ap- 
proximates far-infrared or sub-mm continuum emission from cool 
dust. The green band of the long-wavelength images shows a simple 
approximation to the mid-infrared emission from poly cyclic aro- 
matic hydrocarbons (PAH). The PAHs are assumed to reprocess a 
fixed fraction of the local far-ultraviolet radiation field. Detailed 
calculations (Bakes, Tielens, & Bauschlicher 2001) show that this 
assumption is correct to within a factor of about two for most PAH 
emission features over a broad range of conditions. The local PAH 
density is assumed to be proportional to the neutral plus molecu- 
lar hydrogen density, since there is strong evidence that PAHs are 
destroyed in ionised gas (Bregman et al. 1989; Lebouteiller et al. 
2007). No attempt is made to discriminate between different PAH 
charge states. The blue band of the long- wavelength images shows 
the 6 cm radio continuum emission due to bremsstrahlung in the 
ionised gas. Apart from a slightly different temperature dependence 



this is very similar to the Ha emission, except that it does not suffer 
any dust absorption. 

Starting with the O-star case, we see that the morphological 
appearance shows the typical characteristics of H II regions in turbu- 
lent media found in earlier work (Mellema et al. 2006a; Dale et al. 
2007; Mac Low et al. 2007), namely a fairly irregular structure with 
fingers or pillars, as well as bar-like features at the edge of the H II 
region. Note how the most prominent finger in the bottom edge of 
the H II region does not exactly point towards the source of ionising 
photons. In appearance the simulated Hll regions look strikingly 
like observed H II regions. 

Comparing the magnetic with the non-magnetic case, one sees 
an overall agreement in the structures; the presence of a magnetic 
field does not have a large impact on the global morphology of the 
H II region. This is to be expected as for the O-star case the critical 
time (see Eq. 4) has not been reached by the time most of the volume 
has become ionised. However, the magnetic field does have an effect 
on small scale features in the H II region. In the magnetic case the 
fingers and bars look smoother and broader, due to the magnetic 
pressure, which modifies the radiation-driven implosion of these 
structures (Henney et al. 2009; Ryutov et al. 2005). This behaviour 
is visible in both optical and long-wavelength images, but most 
clearly in the latter. In section 4 we discuss the behaviour of the 
globules in greater depth. 

For the B-star case, the critical time does occur within the 
simulation. However, as is apparent from Fig. 8 this does not lead 
to a significant deformation of the H II region. The reason is that 
magnetic field itself also has a turbulent structure, and there is no 
large-scale field which can impress its direction on the shape of the 
Hll region, as in the test case presented above (see § 2.2.2). Al- 
though the magnetic and non-magnetic results do differ, the overall 
impression is that they are still quite similar. In fact, the difference 
between the O-star and the B-star cases is much larger than between 
the magnetic and non-magnetic cases. For the B-star case, no pil- 
lars are found and the ionised region appears more spherical. The 
edge of the H II region is often more fuzzy, although some sharper 
bar-like features are present. The long wavelength images show the 
reason for this: the H II region is embedded in a thick PDR region 
(with an extent of ~ 30% of the radius of the H II region) which 
erases much of the small-scale density fluctuations present in the 
original cold molecular medium. Close inspection reveals that the 
high-density regions already inside the PDR region photo-evaporate 
and thus lose their large density contrast. In the O-star case these are 
the features which develop into pillars. The fact that the temperature 
in the PDR is closer to 100 K than to 10 K in combination with the 
low ionising flux of the B-star should generally be less conducive to 
the formation of pillars, as shown by Gritschneder et al. (2010). 

Careful examination of the images shows that otherwise the 
effect of the magnetic field is quite similar to that in the O-star 
case: small-scale structures are smoothed out in the magnetised 
case. However, the effect is quite marginal and does not give a clear 
observable diagnostic for the presence of magnetic fields or not. 

3.2.2 Global Properties of the H II Region Evolution 

The initial mean density in the computational box is n = 1000 cm -3 
but the material is distributed inhomogeneously with the densest 
clumps having densities > 10 5 cm -3 . In Figs. 9 and 10 we can 
see how the mean densities in the ionised, neutral and molecular 
components evolve with time and the growth of the H II regions 
around the O and B stars, respectively. In these figures, lines with 
symbols are for the MHD simulations and lines without symbols 
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Figure 7. Simulated optical (top) and long-wavelength (bottom) emission for an H II region around an O star after 200,000 yrs of evolution for MHD (left) and 
purely hydrodynamic (right) simulations. The viewing angle is 6 = 350°, = 350°, where 6 is the polar angle, measured from the x-axis, and is the azimuthal 
angle, measured around the x-axis. Top: Synthetic narrowband optical emission-line images in the light of [Nil] 6584 A(red), Ha 6563 A(green), and O III 
5007 A(blue). Bottom: Synthetic images in the light of 6 cm radio free-free emission (blue), generic PAH (green) and molecular gas column density (red) 



are the purely hydrodynamic case. We can see immediately that the 
magnetic field has a negligible effect on the densities of the different 
components in both the O and B star cases. As the evolution of the 
H II region proceeds, the mean density of the molecular gas in the O 
star simulation increases slowly with time because the lower density 
molecular gas becomes ionised or incorporated into the neutral PDR 



and also because the dense clumps and filaments at the edge of 
the H II region are compressed due to radiation-driven implosions, 
thereby increasing their density. By the end of the simulation, only 
these densest clumps remain, since the majority of the computational 
box has been ionised. In the B star case, the mean density of the 
molecular gas remains roughly constant throughout the simulation. 
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Figure 8. Simulated optical (top) and long-wavelength (bottom) emission for an Hll region around a B star after 10 6 yrs of evolution for MHD (left) and purely 
hydrodynamic (right) simulations. The angles 6 and are defined in the caption to Fig. 7. Top: Synthetic narrowband optical emission-line images in the light of 
[Nil] (blue), Ha (green), and S II (red). Bottom: Synthetic images in the light of 6 cm radio free-free emission (blue), generic PAH (green) and molecular gas 
column density (red) 



This is because the H II region does not break out into regions of 
lower density and remains confined within the molecular gas for the 
duration of the simulation. Also, there is less fragmentation in the 
B star simulation because density inhomogeneities in the neutral 
region of the PDR are smoothed out by photoevaporation flows 



due to FUV radiation before the ionisation front reaches them (see 
Fig. 8). 

By the end of the respective simulations, the mean density in 
the ionised gas around the O star is ~ 100 cm -3 , while that around 
the B star is ~ 10 cm -3 even though the spatial extent is much 
smaller. This is because the latter is a much weaker ionising source. 
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Figure 9. Mean densities of the ionised (solid line), neutral (dashed line) and 
molecular (dotted line) components for the evolving H II region around the 
O star. Lines with symbols are for the MHD simulation, while lines without 
symbols are for the purely hydrodynamic case. 
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Figure 10. Same as Fig. 9 but for the B star. 



In both cases, the density of the neutral PDR gas tends to a roughly 
constant value of 10 3 cm -3 . 

In Figs. 11 and 12 we compare how the ionised, neutral and 
molecular gas component fractions vary throughout the evolution 
of the Hll regions. We calculate both volume and mass fractions 
for each component. The ionised gas in both the O and B star 
cases is thermally dominated and there is no discernible distinction 
between the MHD and purely hydrodynamic results. In the O star 
case, the ionised fraction comes to occupy 90% of the volume but 
only 10% of the mass of the computational box by the end of the 
simulation. In the B star case, the H II region does not manage to 
globally break out of the computational box and so even at the end 
of the simulation the volume occupied by ionised gas is less than 
30%. The mass fraction is < 1% in this case, since the density 
in the ionised gas is low. The neutral gas is most affected by the 
magnetic field for both O and B star simulations. Figs. 11 and 12 
show that the MHD simulation has both higher mass and higher 
volume fractions of neutral gas than the purely hydrodynamic results. 
The neutral component is more important in the B star case than in 
the O star case, reaching ~ 70% of the mass and volume fractions 
after 500,000 yrs and remaining roughly constant thereafter. In the 
O star case, however, the mass fraction remains roughly constant 
for the second half of the simulation while the volume fraction 
decreases sharply. This difference in behaviour can be attributed to 
the greater fragmentation in the O star case (both MHD and HD) 
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Figure 11. Ionised, neutral and molecular gas fractions by volume (solid 
lines) and mass (dashed lines) for the evolving H II region around the O 
star. Lines with symbols are for the MHD simulation, while lines without 
symbols are for the purely hydrodynamic case. 



compared to the B star. The neutral material in the B star simulation 
is distributed in a quite broad, smooth, almost uniform density region 
around the H II region, whereas in the O star simulation the neutral 
region around the photoionised gas is relatively thin but there is 
also neutral material inside the dense clumps and filaments that are 
formed by radiation-driven implosion, which have high density but 
low volume. 

The expansion of an H II region in a clumpy medium is not so 
easy to characterise. In Figs. 13 and 14 we plot a variety of different 
measures of the radius of the photoionised region as functions of 
time, together with the analytical solution obtained using Equation 1 
for a medium of uniform density n = 10 3 cm -3 . We calculate the 
mean radius 2 and also the minimum and maximum radii of the 
ionisation front. For the H II region around the O star, the mean 

2 We define the mean radius as the radius of a sphere with the same volume 
as the real H II region. 
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Figure 12. Same as Fig. 11 but for the B star. 
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Figure 13. Expansion of the H II region with time for the O star. The thick 
grey line represents the analytical solution. Also shown are the mean radius 
(solid lines), the maximum radius (dotted lines) and the minimum radius 
(dashed lines) of the ionisation front. Lines with symbols are for the MHD 
simulation, while lines without symbols are for the purely hydrodynamic 
case. 
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Figure 14. Same as Fig. 13 but for the B star. 



radius closely follows the analytical solution for a long period of the 
simulation, possibly because the filling factor of the dense clumps 
and filaments which retard the ionisation front is small. The maxi- 
mum radius for the ionisation front for the O star simulation quickly 
becomes equal to the distance to the edge of the computational box 
and therefore has no physical meaning after this point. Until it leaves 
the box, the maximum radius represents the direction in which the 
ionisation front is able to expand most quickly, i.e. a path of lowest 
density radially outward from the star. 

For the B star simulation, the mean radius is always below the 
analytical solution. Since this is true for both the MHD and purely 
hydrodynamic simulations, it cannot be due to the magnetic field. 
The reason is that for the B star, the initial clump in which the star 
is embedded has a greater effect on the subsequent evolution of the 
H II region than in the O star case. In fact, we see that the maximum 
radius of the ionisation front in the B star H II region does follow 
the analytical solution. The maximum radius follows the expansion 
of ionisation front in the direction where it was first able to break 
out of the dense clump in which the star was embedded. We have 
already seen that the neutral medium around the H II region in the B 
star case takes on a smooth and homogeneous appearance and that 
the mean molecular and neutral gas densities are roughly constant 
in time. Once the ionisation front has managed to break out of the 
clump, therefore, it follows the expansion law for a uniform medium 
in this mean density. In the majority of directions, however, the 
ionisation front takes much longer to break out of the initial clump 
and so the mean ionisation front radius is retarded compared to the 
analytical solution. 

In both O and B star simulations, the minimum radius of the 
ionisation front indicates the presence of fingers and clumps of neu- 
tral material within the H II region. For the O star simulation we 
have already mentioned how the magnetic field suppresses fragmen- 
tation, and this is reflected in the graph of minimum radius with time: 
fingers and globules formed in the purely hydrodynamic simulation 
are denser and survive longer closer to the star compared to the 
MHD simulation. In the B star case there is much less fragmentation 
and we see no fingers and clumps of neutral gas within the H II 
region and the minimum radii of both MHD and hydrodynamic 
simulations seen in Fig. 14 bear this out. 

In Figs. 15 and 16 we show the mean radial velocities for 
the three gas components. Again, we compare the results of both 
MHD (lines with symbols) and purely hydrodynamic simulations. 
Beginning with the ionised gas component, we see that for the O star 
the rms velocities reach a peak of about 8 km s _1 , and then decreases 
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Figure 15. Mean gas radial velocities within the ionised, neutral and molec- 
ular gas components for the O star simulation. The thick grey line represents 
the analytical mean radial velocity. Also shown are the rms radial velocities 
(solid lines) and the mean radial velocity (dashed lines) of the gas. Lines 
with symbols are for the MHD simulation, while lines without symbols are 
for the purely hydrodynamic case. 



very slowly, which is what we saw in our previous paper (Mellema 
et al. 2006a). The mean radial velocity, on the other hand, peaks at 
8-10 km s _1 then declines quite quickly. The rms velocities are due 
to interacting photoevaporated flows within the H II region, which 
lead to gas flowing inwards as well as outwards (Henney 2003). The 
mean radial velocity traces the general expansion of the photoionised 
gas. We see that for the non-uniform medium, the initial expansive 
motions are more rapid than for the simple analytical model but this 
is just because the gas first begins to expand into regions of lowest 
density (Shu et al. 2002). There is a small difference between the 
rms velocities in the MHD and hydrodynamic cases, with the latter 
being marginally higher. This could be due to the greater degree 





250 500 750 1000 1250 1500 

Time, 1000 yr 

Figure 16. Same as Fig. 15 but for the B star. 



of fragmentation in the hydrodynamic simulations which, in turn, 
leads to stronger photoevaporated flows. 

In the B star simulation the rms velocities are similar to the 
O star case just discussed, which suggests that photoevaporated 
flows are important in this H II region, too. The pure hydrodynamic 
simulation shows slightly higher velocities than the MHD case, 
which indicates that the magnetic field plays some role in reducing 
the importance of photoevaporated flows within the ionised gas. The 
mean radial velocities show the same sort of ringing seen in the test 
problem of H II region expansion in a uniform density medium (see 
Fig. 2). The period of the oscillations is of order 2 X 10 5 yrs, which 
is consistent with the sound crossing time of the photoionised region. 
Note that in both O and B star simulations, the actual mean radial 
velocities are considerably higher than the corresponding analytical 
values for expansion in a uniform medium with density equal to the 
mean density of the initial computational box (n = 1000 cm -3 ). 

The neutral gas velocities show more separation between hy- 
drodynamic and MHD results, with the hydrodynamic velocities 
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being consistently higher in both mean and rms velocities. Inter- 
estingly, the different simulation velocities appear to depart from 
each other at a specific point in time, which is different depending 
on whether the mean or rms velocity is considered. This is true for 
both O and B star simulations. In the O star case, after the initial 
acceleration, the neutral gas velocities reach 4-5 km s _1 , while for 
the B star the velocities are 1-2 km s _1 lower and show evidence for 
ringing as for the ionised gas. 

The initial molecular gas velocities reflect the Mach 10 su- 
personic turbulence of the underlying density distribution. At early 
times, there is evidence of residual infalling motions, i.e. negative 
radial velocities, again from the initial conditions. These persist 
longer in the B star simulations than in the O star simulations. This 
is because the O star rapidly ionises out to the edge of the com- 
putational box, and the molecular material which remains is being 
accelerated radially away from the central star in clumps that are 
experiencing the rocket effect. In the B star case, most of the molec- 
ular material surrounding the H II region remains undisturbed, and 
so the imprint of the initial conditions remains in this gas until the 
end of the simulation. 

In summary, we find slight differences in the expansion of H II 
regions between MHD and purely hydrodynamic simulations but 
the greatest differences are due to the different ionising and FUV 
fluxes of the central O or B star. Most of the differences between 
the MHD and hydrodynamic cases are due to the lesser amount of 
fragmentation in the former, indicating that magnetic fields provide 
some support to the neutral gas against radiation-driven implosions. 

3.3 Magnetic Quantities 

Even though we find that the presence of a magnetic field only has 
small effects on the growth of the H II region, this does not mean 
that the field is uninteresting or unimportant. For one thing, the 
field is important in the still neutral medium, where it constitutes a 
significant fraction of the total energy. 

In Fig. 17 we show the magnetic field and the distribution of 
the photoionised gas for the B star case in the same format as Figs. 3 
and 4. We can see that the magnetic field looks very different. The 
initial conditions were the result of an MHD turbulence calculation 
and in the molecular gas (coloured grey/white in the lower panels of 
Fig. 17) we can see that the field appears to be randomly oriented. 

In the warm, neutral (PDR) gas (coloured purple in the lower 
panels of the figure) there is strong evidence that the magnetic 
field is being aligned parallel to the ionisation front, particularly 
in the lower left region of the x-y and x-z plots. At some points in 
this image, the ionisation front has left the grid, and the periodic 
boundary condition in this simulation means that ionised material 
exits one side of the grid and reappears at the opposite boundary, 
where it will usually recombine because of the high column density 
through the molecular and neutral gas to the star at this position. At 
this point in the simulation, this is not affecting the global evolution 
of the H II region and surrounding warm, neutral gas. 

The photoionised region has a very low magnetic field intensity, 
showing that the field has been pushed aside by the expansion of the 
H II region. Photoevaporated flows from density inhomogeneities 
in both the neutral and molecular gas drag the magnetic field with 
them into the photoionised region in long filamentary structures. 
However, there is little mass associated with these flows. In the x-y 
and x-z central planes there is a large region of partially ionised 
gas. This is also in the direction where the ionisation front has left 
the grid and could be due to a density or velocity gradient, which 
results in the ionisation front being unable to keep up with the gas, 



since it is no longer preceded by a neutral shock in this part of the 
computational grid. 

In order to analyze the role of the magnetic field in our three 
components, ionised, neutral and molecular, we present colour plots 
showing the joint distribution of various quantities in these three 
components. Since the B-star and O-star results are quite compara- 
ble, we only show the B-star results. 

We start by considering joint distribution of the magnetic field 
strength against the density. Fig. 18 shows this for the case of the 
B-star at t = 10 6 years. The distribution of the ionised material is 
shown in blue and displays a wide range of field strengths (0.1 to 10 
yt/G) across a narrow range of densities (10 to 20 cm -3 ; in the case of 
the O-star these values are some 5 times higher). The diagonal lines 
show the Alfven speeds. If the sound speed in the gas is smaller 
than v A the flow will be magnetically dominated. As can be seen, 
the ionised material does not reach Alfven speeds of 10 km s _1 , 
the typical sound speed in the ionised medium. Consequently the 
ionised flow is never magnetically dominated, consistent with the 
results found in the previous sections. 

Another way to represent the role of the magnetic field is to 
plot the joint distribution of magnetic and ram (turbulent) pressure, 
as is done in Fig. 19. In full equipartion, the turbulent, magnetic 
and thermal pressures should all be equal. Since the H II region is 
expanding through its thermal pressure, it is the thermal pressure 
which is the overall dominating one. Fig. 19 shows that of the 
other two pressure components, turbulent pressure mostly dominates 
over magnetic pressure. In only a few places in the H II region the 
magnetic pressure dominates. 

In the neutral and molecular material (represented by green 
and red in the joint distribution figures) the situation is different. 
The joint distribution of density and magnetic field strength, Fig. 18, 
shows a wide range of densities (10 2 -10 4 cm -3 ) and a narrower 
range of magnetic field strengths (mostly around 10-30 //G, although 
some areas, principally neutral ones, have magnetic fields as weak as 
1 iiG). The molecular material generally has higher densities than the 
neutral material. Since the sound speed in the neutral and molecular 
regions is 1 km s _1 or less, they are mostly magnetically rather 
than thermally dominated. In the distribution of pressures, Fig. 19, 
one sees that the neutral and molecular material clusters around 
equipartition of the magnetic and turbulent pressures. The molecular 
material has regions where the turbulent pressure dominates, and 
others where the magnetic pressure dominates. The neutral (PDR) 
material has generally somewhat lower pressures and is closer to 
equipartition. 

Fig. 18 can be compared to observationally derived values 
for the magnetic field and densities in H II regions and their sur- 
roundings. Results of Harvey-Smith et al. (in preparation) on large 
Galactic H II regions show ranges in magnetic field strengths and 
typical densities for the ionised, neutral and molecular components 
very comparable to what we find in our simulation results. 



4 DISCUSSION 

4.1 Comparison with observations: RCW 120 

The H II region RCW 120 and its surrounding medium have recently 
been extensively studied at near- to far-infrared wavelengths in the 
context of triggered star formation (Zavagno et al. 2007; Deharveng 
et al. 2009; Anderson et al. 2010; Martins et al. 2010; Zavagno et al. 
2010). This small (3.5 pc diameter at 1.35 kpc distance, Zavagno 
et al. 2007) H II region appears to have a single ionising source, 
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t = 1.00 Myr 



Figure 17. Magnetic field and photoionised gas distribution in the central x-y (left), x-z (centre) and y-z planes for the B star simulation after 10 6 yrs. Colours 
have the same meaning as in Fig. 3. 



identified from VLT-SINFONI near-IR spectroscopy and spectral 
line fitting to stellar atmosphere models as an O6-O8V/III star 
(Martins et al. 2010) with an effective temperature of T eS = 37.5 + 
2 kK and ionising photon rate of log <2h = 48.58 + 0.22. These 
observationally derived parameters are remarkably similar to those 
we have adopted for our O star simulations and so an opportunity 
for a direct comparison of RCW 120 with our results presents itself. 

Morphologically, the observations and simulation look very 
similar. In particular, our Fig. 7, corresponding to the O star simula- 
tion after 200,000 yrs of evolution, is the same size as the observed 
bubble RCW 120. The upper panels of the figure show the optical 
emission (with extinction from dust), while the lower panels show 
synthetic 6 cm radio free-free emission from the ionised gas (blue), 
generic PAH emission from the PDR (green) and cold molecular gas 
column density (red). Figs. 1 and 2 of Deharveng et al. (2009) show 
the Ha emission from the ionised gas, the 8 jim emission from PAHs 
in the PDR, 24 /im emission from warm dust in the photoionised re- 
gion, and 870 fim cold dust emission from the neutral and molecular 
gas. If the reader mentally rotates the simulated images so that the 
lower left corner moves to the top centre, then a certain correspon- 
dence between simulations and observations can be imagined. The 
ionised gas emission fills the central cavity. In the optical, the images 
appear criss-crossed by dust lanes. Some of these are foreground, 
others are the result of projection of filaments and ridges deeper 
inside the H II region. The PAH emission comes from a thin shell 



around the periphery of the H II region, with projection seeming 
to put some PAH emission in the interior of the H II region, as can 
be seen in both the simulation and the observations. The simulated 
PAH emission presents the same sort of irregularities, filaments and 
clumps as are seen in the observations. 

The initial total mass in the computational box is ~ 2000 M Q . 
After 200,000 yrs of evolution, the mass distribution, as seen in 
Fig. 11, is approximately 5% ionised gas, 55% neutral gas, and 40% 
molecular gas. Some of the molecular gas is undisturbed material 
at the edge of the computational box but quite a large proportion 
must be in the clumps and filaments in the H II region and PDR, 
where it becomes compressed by radiation-driven implosion. The 
H II region has influenced virtually all of the computational domain 
by this time, and the ~ 2000 M Q of neutral and molecular material is 
comparable to the 1 100-2100 M Q derived by Zavagno et al. (2007) 
and Deharveng et al. (2009) from the 870 jim cold dust emission. 
This dust will be present in both the neutral and molecular gas. 
Hence, we can surmise that the average density in the vicinity of 
RCW 120 is similar to that in our computational box, that is (uq) = 
10 3 cm -3 , though, of course, the density distribution is far from 
homogeneous. From this, we could go so far as to postulate that 
the age of RCW 120 must be similar to the time depicted in our 
simulation, that is, of order 200,000 yrs. Martins et al. (2010) were 
unable to assign an age to RCW 120 except to say that it must be 
younger than 5 Myr. 
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Bstar-ep, t = 1.000 Myr 
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Figure 18. Magnitude of B-field (Gauss) against number density (cm -3 ) for 
volume-weighted joint distributions in the case of the B-star at t = 10 6 years. 
The diagonal lines correspond to constant Alfven speeds of 10 and 1 km s -1 . 
The three colours represent ionised (blue), neutral (green) and molecular 
(red) material. 
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Figure 19. Magnetic pressure versus turbulent pressure for volume-weighted 
joint distributions in the case of the B-star at t = 10 6 years. Note that the 
scales are the same on both x and y axes. The diagonal line shows equal 
pressures. The three colours represent ionised (blue), neutral (green) and 
molecular (red) material. 



Although radiation-driven implosion does enhance the density 
of photoionised globules at the periphery of the simulated Hll 
region, the fact that we do not include self gravity in our simulations 
means that we are unable to model triggered star formation in the 
neutral shell around the H II region. However, it is quite likely that 
radiation-driven implosion could trigger gravitational collapse in 
clumps that are already on the point of forming stars. 

There are other potentially important physical processes that we 
do not include in our models. Dust grains are an important compo- 
nent of the interstellar medium in star-forming regions. They absorb 
ionising photons and both observations (e.g., Wood & Churchwell 
1989) and theory (e.g., Arthur et al. 2004) suggest that more than 
50% of all the ionizing photons are absorbed by dust within the H II 
regions in the initial stages of evolution when the ionised density 
is high. Radiation pressure on dust grains can form a central cavity 
in the Hll region (Mathews 1967; Inoue 2002), which can be as 
large as 30% of the ionisation front radius. Such cavities have been 
known for a long time from optical observations of H II regions 
(e.g., Menon 1962) and are also seen in more recent infrared ob- 
servations (e.g., Watson et al. 2008). Krumholz & Matzner (2009) 
have studied the effect of radiation pressure on the dynamics of H II 
regions and conclude that it is only an important effect for massive 
star clusters and not for H II regions around individual massive stars. 
Stellar winds are also to be expected from stars whose effective 
temperatures are greater than 25,000 K and these winds will provide 
an alternative mechanism for evacuating a central cavity in the dust 
and gas distribution in the H II region. The mechanical luminosity of 
the stellar wind is converted into thermal pressure by shock waves 
(Dyson & Williams 1997) and this could affect the dynamics of the 
H II region. 

4.2 Predicted maps of projected magnetic field 

Figs. 20 to 23 show visualisations of the projected integrated B-field 
(line-of-sight and plane-of-sky components) from our simulations. 
These visualisations represent generic idealised versions of maps 
that can be obtained from various observational techniques. The 
line-of-sight field component can be determined from the Faraday 
rotation measure (Heiles et al. 1981) for ionised gas or from Zee- 
man spectroscopy (Brogan & Troland 2001) for neutral/molecular 
gas. The plane-of-sky field components can be determined from 
observations of polarized emission or absorption (Kusakabe et al. 
2008) from aligned spinning dust grains. In all of these techniques, 
the determination of B is not straightforward and relies on many 
auxiliary assumptions. In particular, it is very difficult to determine 
the absolute magnitude of the plane-of-sky field components except 
via statistical techniques such as the Chandrasekhar-Fermi method 
(Chandrasekhar & Fermi 1953; Heitsch et al. 2001; Ostriker et al. 
2001; Falceta-Goncalves, Lazarian, & Kowal 2008). Furthermore, 
many of the techniques do not uniformly sample all regions along 
the line of sight, but are biased towards regions with particular phys- 
ical conditions. We therefore caution against direct comparison of 
observational results with our maps in any but a qualitative sense. 
The maps are nonetheless useful in showing an overview of the 
magnetic field geometry in our simulations. In the following discus- 
sion, we describe locations on the projected map using the points of 
the compass, assuming that north is up (positive y) and east is left 
(negative x). 

The most striking aspect of Figs. 20 and 21 is the large-scale 
order that is apparent in the projected magnetic field. This is particu- 
larly visible in the neutral- weighted maps (top-right panels), where 
it can be seen that the field is frequently oriented parallel to the 
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Figure 20. Projected magnetic field of O star model at a time of t = 200, 000 years, showing the full simulation cube (extent of 4 pc along each axis). The 
horizontal and vertical graph axes correspond to the simulation x and y axes, respectively. The three grayscale images (top-left, top-right, and bottom-right 
panels) show the line-of-sight integral of the line-of-sight component of the magnetic field, weighted by the ionised (top-left), neutral (top-right), or molecular 
(bottom-right) gas density. The superimposed red vectors indicate the magnitude and direction of the line-of-sight integral of the plane-of-sky components of the 
magnetic field, using the same weightings. The integrals of the plane-of-sky components were carried out in a Stokes QU frame, as described in the appendix. 
Blue contour lines show the column density of molecular gas, to aid comparison between the panels. The bottom-left panel shows a composite image of the 
column densities of molecular (red), neutral(green), and ionised (blue) gas. The normalisation varies between the panels. For the column densities (bottom-left 
panel and contours) the maximum plotted values in units of H nuclei per cm 2 are 3.7 x 10 21 (ionised), 5.1 x 10 22 (neutral), and 1.3 x 10 23 (molecular). For the 
projected magnetic fields, the maximum plotted values in units of //G times H nuclei per cm 2 are 5.2 x 10 22 (ionised, top left), 2.1 x 10 24 (neutral, top right), 
1.7 x 10 24 (molecular, bottom right). 
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Figure 21. Projected magnetic field of B star model at a time of t = 10 6 years, showing the full simulation cube (extent of 4 pc along each axis). Panels as in 
Fig. 20. For the column densities (bottom-left panel and contours) the maximum plotted values in units of H nuclei per cm 2 are 2.1 x 10 20 (ionised), 4.6 x 10 22 
(neutral), and 3.4 x 10 22 (molecular). For the projected magnetic fields, the maximum plotted values in units of /uG times H nuclei per cm 2 are 1.0 x 10 21 
(ionised, top left), 7.5 x 10 23 (neutral, top right), 7.2 x 10 23 (molecular, bottom right). 



large-scale ionisation front, forming a ring around the Hll region. 
The effect is particularly strong along the north and south borders 
of the region because of the net positive field along the x-axis (see 
§ 3.1), but it is also seen to a lesser extent along the east and west 
borders, despite the fact that the mean [/-component of the field is 
zero. This is because of the fast-mode MHD shock that is driven 
into the surrounding gas by the expanding H II region and which 
compresses both the gas and the field, tending to bend the field lines 
so that they are more closely parallel to the shock than they were 
in the undisturbed medium. A large-scale pattern is much harder 



to discern in the molecular-weighted maps (bottom-right panels), 
partly because the molecular column density is much less smoothly 
distributed, being concentrated in globules and filaments. 

In most dense filaments, the field direction is parallel to the 
long axis of the filament, as can be seen particularly clearly in 
Fig. 22, which shows a detail view of the dense photoevaporating 
globule found to the south in Fig. 20 and which is fed by multiple 
neutral/molecular filaments. In the molecular gas at the head of the 
globule, the magnetic field is bent into a hairpin shape. The B-field 
in the ionised gas at the head of the globule tends to be oriented 
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Figure 22. Projected magnetic field of O star model at a time of t = 200, 000 years, showing a 4x zoom centred on the southern dense globule (extent of 1 pc 
along each axis). Panels as in Fig. 20. For the column densities (bottom-left panel and contours) the maximum plotted values in units of H nuclei per cm 2 are 
3.7 x 10 21 (ionised), 2.7 x 10 22 (neutral), and 1.3 x 10 23 (molecular). For the projected magnetic fields, the maximum plotted values in units of /uG times H 
nuclei per cm 2 are 2.6 x 10 22 (ionised, top left), 1.0 x 10 24 (neutral, top right), 2.2 x 10 24 (molecular, bottom right). 



perpendicular to the ionisation front, as is the case with almost all 
the dense globules visible in Fig. 20. The same is true along much 
of the bar-like feature to the west of the globule. On the other hand, 
in a few regions, such as along the filament that extends south from 
the globule, the B -field in the ionised gas lies along the ionisation 
front. 

Fig. 23 shows a detail view of the same dense southern globule, 
but from the B star simulation shown in Fig. 21. The molecular 
gas shows a similar field pattern to that seen in the O star case: the 
field goes up one of the feeding filaments and down the other, with 



a hairpin bend in between. Although this can be seen clearly in 
the molecular gas, the neutral and ionised gas show very different 
patterns, with magnetic field vectors that are generally perpendicular 
to the filament and continuous with the large-scale field pattern in 
the region. This is because the total ionised and neutral columns are 
dominated by diffuse material along the line of sight, rather than by 
material associated with the globule and filament. 
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Figure 23. Projected magnetic field of B star model at a time of t = 10 6 years, showing a 4x zoom centred on the southern dense globule (extent of 1 pc 
along each axis). Panels as in Fig. 20. For the column densities (bottom-left panel and contours) the maximum plotted values in units of H nuclei per cm 2 are 
1.7 x 10 20 (ionised), 1.5 x 10 22 (neutral), and 3.2 x 10 22 (molecular). For the projected magnetic fields, the maximum plotted values in units of /uG times H 
nuclei per cm 2 are 5.1 x 10 20 (ionised, top left), 3.8 x 10 23 (neutral, top right), 3.6 x 10 23 (molecular, bottom right). 



4.3 Effects of magnetic field on globule formation and 
evolution 

It is interesting to compare the properties of the globules generated 
in our turbulent simulations with the results of previous detailed 
studies of the photoionisation of isolated dense globules (Williams 
2007b; Henney et al. 2009; Mackey & Lim 2010b). The principal 
findings of the earlier studies are that a sufficiently strong magnetic 
field (J3 < 0.01 in the initial neutral globule) will produce important 
qualitative changes in the photoevaporation process. Depending 



on the initial field orientation with respect to the direction of the 
ionising photons, either extreme flattening of the globule may occur 
(perpendicular orientations) or the radiative implosion of the globule 
may be prevented (parallel orientations). For all orientations the 
ionised photoevaporation flow cannot freely escape from the globule, 
leading to recombination at late times. On the other hand, a weaker 
field (J3 0.01 in the initial globule) leads only to a moderate 
flattening of the globule and does not prevent the free escape of the 
ionised photoevaporation flow. 
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In our turbulent simulations, the mean initial value of J3 is 
0.032 (§ 3.1), which is intermediate between the two cases dis- 
cussed above and might lead one to suspect that magnetic effects on 
the evolution of globules should be substantial. However, this is not 
the case. A careful examination of the three-dimensional globule 
morphologies in our simulations shows no evidence of magnetically 
induced flattening. Although many globules do show asymmetries, 
this is true equally of our non-magnetic simulations and is pre- 
sumably due to their irregular initial shape and internal turbulent 
motions. The only difference in the globule properties between our 
non-magnetic and magnetic simulations is that the very smallest 
globules do not seem to form in the magnetic case. 

In order to explain this apparent discrepancy with earlier work, 
it is necessary to examine in more detail the distribution of magnetic 
field in the initial conditions of our turbulent simulations. It turns out 
that the dense filaments of molecular gas (from which the globules 
will ultimately form) tend to be much less magnetically dominated 
than the more diffuse gas, typically showing J3 > 0.1. Furthermore, 
these filaments already show a magnetic geometry similar to that 
described above (§ 4.2), with the field running along the long axis 
of the filament and with the field changing sign as one moves across 
the short axis. This is very different from the initial conditions 
assumed in the earlier globule photoevaporation studies, which 
were a uniform magnetic field that threaded a spherical globule 
(e.g., Fig. 1 of Henney et al. 2009). As a result, the compressed 
globule heads in our turbulent simulations show approximately 
equal thermal and magnetic pressures (J3 ~ 1) and the magnetic 
effects are very modest. One caveat to this result is that numerical 
diffusion due to our limited spatial resolution may be producing non- 
physical magnetic reconnection in the globule head when oppositely 
directed field lines are forced together. Higher resolution studies 
of globule implosion with realistic initial field configurations are 
required in order to clarify this and these should also include self 
gravity, which may be important during the phase of maximum 
compression (Esquivel & Raga 2007). 



5 CONCLUSIONS 

We have performed radiation-magnetohydrodynamic simulations of 
the formation and expansion of H II regions and PDRs around an 
O star and a B star in a turbulent magnetised molecular cloud on 
scales of up to 4 parsec. Our principal conclusions are as follows: 

(i) The expansion of the H II region is little affected by the pres- 
ence of the magnetic field, since the thermal pressure of the ionised 
gas dominates the dynamics on the timescales of our simulations 
(§ 3.2.2). 

(ii) The O star simulations produce greater fragmentation and 
denser clumps and filaments around the periphery of the H II region 
than the B star case (§ 3.2.1, Figs. 7 and 8). 

(iii) For B stars the non-ionising far ultraviolet radiation plays 
an important role in determining the morphology of the region. 
H II regions around such stars are surrounded by a thick, relatively 
smooth shells of neutral material (PDR), ~ 30% of the bubble radius 
(e.g., lower panels of Fig. 17). In the O star simulations, the PDR is 
thinner and more irregular in shape (e.g., Fig. 7). 

(iv) The resemblance at optical and longer wavelengths of our 
simulations to observed bubbles is striking (§ 4.1 and Deharveng 
et al. 2009, 2010). Our ~ 2 pc radius bubbles are a typical size 
compared to bubbles at known distances in the GLIMPSE surveys 
(Churchwell et al. 2006, 2007), and so comparisons are meaningful. 



(v) The expanding H II region and PDR tend to erase pre-existing 
small-scale disordered structure in the magnetic field, producing 
a large-scale ordered field in the neutral shell, with orientation 
approximately parallel to the ionisation front (top panels of Fig. 17 
and top-right panels of Figs. 20 and 21). 

(vi) Dense evaporating globules, pillars, and elephant trunk struc- 
tures tend to be fed by two or more neutral/molecular filaments, with 
magnetic fields running along their length (right panels of Fig. 22 
and bottom-right panel of Fig. 23). The field geometry in the neutral 
and molecular gas at the bright head of the structure tends to be of a 
hairpin shape. 

(vii) The weak magnetic field in the ionised gas also shows an 
ordered structure. On the largest scales and at the latest times, it 
tends to align (albeit weakly in the O star case), with the mean 
field direction of the simulation's initial conditions (Fig. 21, top- 
left panel). On smaller scales, there is a general (but not universal) 
tendency for the field in ionised gas to be oriented perpendicular to 
the local ionisation front. This tendency is more pronounced in the 
O star simulation (Fig. 20, top-left panel) and in globule evaporation 
flows for both simulations (top-left panels of Figs. 22 and 23). 

(viii) The highest pressure compressed neutral/molecular gas 
shows approximate equipartition between thermal, turbulent, and 
magnetic energy density, whereas lower pressure gas (either neutral 
or molecular) tends to separate into, on the one hand, magnetically 
dominated, quiescent regions, and, on the other hand, demagnetised, 
highly turbulent regions (§ 3.3, Fig. 19). The lower pressure gas also 
separates into low-/?, magnetically dominated regions (which are 
largely molecular) and high-/?, thermally dominated regions (which 
are largely neutral). The ionised gas, on the other hand, always 
shows approximate equipartition between thermal and turbulent 
energies, but with the magnetic energy being lower by 1 to 3 orders 
of magnitude. 

(ix) Velocity dispersions in the ionised gas of 7-9 km s _1 are 
maintained for the entire duration of all our simulations (§ 3.2.2, 
Fig. 15 and 16). This is 5 to 10 times higher than the value that 
would be predicted by expansion in a uniform medium. At early 
times (t < 100, 000 yr for the O star, or t < 300, 000 yr for the B star), 
this dispersion is mainly due to radial champagne-flow expansion as 
the H II region escapes from its natal clump. At later times, the net 
radial expansion of ionised gas subsides, but the velocity dispersion 
is maintained by inwardly-directed photoevaporation flows from 
globules and pillars. 
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APPENDIX A: INTEGRATION OF PLANE-OF-SKY 
FIELD COMPONENTS 

Polarisation-based methods for measuring the plane-of-sky com- 
ponents of the magnetic field cannot distinguish the sense of the 
field and are degenerate between a position angle 6 and a position 
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angle 6 + 180°. It is therefore not sufficient to simply integrate B x 
and B y along the line of sight with appropriate weighting. Neither 
is it sufficient to integrate B 2 X and B 2 , since this will always give 
a result in the first quadrant (B x , B y > 0). In order to calculate 
our projected magnetic field maps (§ 4.2), we therefore adopt a 
"Stokes parameter" approach (Chandrasekhar 1960), whereby the 
plane-of-sky components B x (x,y,z) and B y (x,y,z) are first trans- 
formed to (2, £/)-space as B Q = B cos 26, B v = 5 sin 20, where 
B = (B 2 X + B 2 y ) 1/2 and 6 = tan" 1 (B y /B x ). Then the line-of-sight in- 
tegration is performed as M QfU (x, y) = f Bq^(x, y, z) n dz, where n 
is the relevant density (ionised, neutral, or molecular). Finally, the 
components of the projected field are transformed back to (x, y)- 
space: M x = Mcos#, M y = Msin#, where M = {M 2 Q + M^) 1/2 and 
6 = 0.5 tan" 1 (M v /Mq). 
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